%Simulate the future state of each hospital by starting from the actual
%state each year, including variable cost savings from the dominant local
%vendor.

tic

choice=1:13;
ns=numel(choice);
beta=0.95;
tbar=3;
%The alpha and gamma are the coefficients from the policy function for
%endogenous hospitals (stand-alone hospitals) that haven't adopted EMR.
%alpha=[0.0163	0.665	0	0	0.476]; %the coefficients from the policy function (vary across choices and hospitals): vmkt vmsh vlag vsys vmkt*big

alpha=[0.176	0.930	1.117	1.205 0 0 0 0  0.0000195];  %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q vdom*yrsince

mktcat_alpha=[1.010	1.964	4.369];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_alpha=[zeros(numel(mktcat_alpha),1),repmat(mktcat_alpha',1,ns-1)];

gamma2=[0.00131	-0.662	-15.54	0.321	0.0959	-7.586]; %the following gamma's are coefficients for the choice-fixed-hospital-varied variables (from estimating the policy function based on the entire sample)
gamma3=[0.00235	0.323	0.875	2.571	0.532	-12.58];  %the order is as follows: preadp(indicate whether adopting or not) bdtot nprofit teaching HHI totpop65 constant 
gamma4=[-0.0107	0.0575	-15.25	2.081	0.224	-6.328];
gamma5=[-0.0113	0.0192	-15.23	1.046	-0.113	-2.999];
gamma6=[0.00874	1.243	-16.39	-1.594	0.0761	-8.910];
gamma7=[0.00283	0.968	1.867	0.0977	-0.00211	-7.356];
gamma8=[0.00925	1.515	-0.894	-0.747	-0.345	-5.642];
gamma9=[-0.00585	-0.396	-15.24	0.383	0.0320	-4.986];
gamma10=[0.00515	0.283	-1.816	2.660	0.271	-9.112];
gamma11=[0.00494	1.089	-1.437	1.722	0.406	-11.26];
gamma12=[0.000429	0.488	-1.904	2.269	0.310	-8.071];
gamma13=[0.00131	0.501	0.0869	2.505	0.293	-8.797];


Gamma=[zeros(size(gamma2))' gamma2' gamma3' gamma4' gamma5' gamma6' gamma7' gamma8' gamma9' gamma10' gamma11' gamma12' gamma13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)


%The coefficients for stand-alone hospitals that have adopted EMR.
tau=[-0.139	0.583	1.180	0.487	6.146	6.179	5.937	5.636	0.0000452]; %the coefficients from the policy function (vary across choices and hospitals): vmkt vmsh vlag vsys vlag*big vdom*yrsince
mktcat_tau=[0	0	0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

zeta2=[0 0 0 0 0 0]; %the following gamma's are coefficients for the choice-fixed-hospital-varied variables (from estimating the policy function based on the entire sample)
zeta3=[0.00108	-1.399	0.219	0.0164	0.0614	0.260];  %the order is as follows: preadp(indicate whether adopting or not) bdtot profit teaching HHI totpop65 constant 
zeta4=[-0.0129	0.867	-10.96	1.740	0.117	0.887];
zeta5=[-0.0194	-1.076	-11.06	0.713	-0.0377	3.168];
zeta6=[0.000848	-1.168	1.119	2.595	0.322	-4.058];
zeta7=[0.00139	-2.752	0.901	1.257	0.233	-1.803];
zeta8=[0.00186	0.208	1.099	4.626	0.853	-12.92];
zeta9=[-0.0107	-1.079	-10.80	-0.931	-0.274	5.019];
zeta10=[0.000679	-2.337	0.246	-0.406	-0.269	4.118];
zeta11=[0.00126	-0.758	0.0698	0.554	0.0997	-1.110];
zeta12=[-0.00109	-2.331	-1.498	1.158	0.172	0.135];
zeta13=[-0.000645	-0.0531	0.575	0.591	0.119	-1.816];

Zeta=[zeros(size(zeta2))' zeta2' zeta3' zeta4' zeta5' zeta6' zeta7' zeta8' zeta9' zeta10' zeta11' zeta12' zeta13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)



%The sigma and eta are the coefficients from the policy function based on
%the sample of exogenous hospitals (affiliated hospitals) that haven't adopted EMR.
sigma=[0 0 0 2.888]; %the coefficient for vsysh 

mktcat_sigma=[0.717	1.513	4.089];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_sigma=[zeros(numel(mktcat_sigma),1),repmat(mktcat_sigma',1,ns-1)];

eta2=[0.00125	-0.612	-0.846	2.580	0.483	-10.92]; %the order is as follows: preadp(indicate whether adopting or not) bdtot nprofit perc_mcr hhi totpop65 constant
eta3=[0.00178	0.212	-0.984	0.381	0.120	-5.473];
eta4=[-0.00891	-0.556	2.063	0.333	0.215	-7.360];
eta5=[-0.0128	-0.253	0.446	1.242	-0.0483	-3.968];
eta6=[0.00121	1.708	-2.342	4.812	0.990	-18.64];
eta7=[0.000986	2.948	-1.246	1.203	0.152	-8.513];
eta8=[0.00206	2.089	-0.481	5.713	0.846	-18.29];
eta9=[-0.00247	-1.280	-2.481	1.469	-0.0275	-2.685];
eta10=[0.00174	0.446	-1.106	1.721	0.328	-8.198];
eta11=[0.000201	0.0966	-3.307	0.200	0.101	-3.988];
eta12=[-0.000114	1.090	-1.144	1.231	0.126	-6.022];
eta13=[0.00135	0.772	-1.957	0.531	-0.0359	-4.213];

Eta=[zeros(size(eta2))' eta2' eta3' eta4' eta5' eta6' eta7' eta8' eta9' eta10' eta11' eta12' eta13'];



%The coefficients for the affiliated hospitals that have adopted EMR.
rho=[0 0 4.342   2.259]; %the coefficient: vmkt vmsh vlag vsys

mktcat_rho=[0    0   0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

lambda2=[0 0 0 0 0 0]; %the order is as follows: preadp(indicate whether adopting or not) bdtot perc_mcr perc_mcd hhi totpop65 constant

lambda3=[0.00181	4.622	3.286	-1.681	-0.167	-0.178];
lambda4=[-0.00498	4.366	3.725	-2.979	-0.648	5.693];
lambda5=[-0.00652	5.133	2.110	-5.853	-1.171	11.33];
lambda6=[0.00237	4.250	1.534	-1.598	0.0581	-3.114];
lambda7=[0.00281	4.065	-0.713	-3.622	-0.377	3.603];
lambda8=[0.00128	2.802	0.0644	-1.662	-0.134	-0.127];
lambda9=[-0.00445	4.364	4.916	-1.729	-0.357	2.017];
lambda10=[0.000729	4.723	5.322	-3.118	-0.348	1.585];
lambda11=[0.0000586	1.672	2.441	-0.909	-0.445	3.386];
lambda12=[-0.000413	4.372	3.812	-1.539	-0.316	1.652];
lambda13=[0.000687	6.875	4.990	-1.436	-0.134	-3.047];
Lambda=[zeros(size(lambda2))' lambda2' lambda3' lambda4' lambda5' lambda6' lambda7' lambda8' lambda9' lambda10' lambda11' lambda12' lambda13'];



bed25=25;
bed50=80;
bed75=193;
%ns=1;       %will take 516 secs if ns=500
T=100;
%s=1;
textFileName3=sprintf('fs%d.csv',s);  %the output file
textFileName4=sprintf('ddom%d.csv',s);   %the output file
textFileName5=sprintf('eerr%d.csv',s);   %the output file
textFileName6=sprintf('mks%d.csv',s);   %the output file

rng(s);
FSTA=[];
DDOM=[];
EERR=[]; %### UPDATED on 03/16/2017: the output including the probability as the expected error
MKTS=[]; %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation


for k=6:9
%k=7;
textFileName1a=sprintf('prdprob_a%d.csv',k); %the input file
textFileName1na=sprintf('prdprob_na%d.csv',k); %the input file
textFileName2a=sprintf('syschs_a%d.csv',k);  %the input file
textFileName2na=sprintf('syschs_na%d.csv',k);  %the input file
J1=dlmread(textFileName1na);  % ### UPDATED on 04/21/2020---the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName1a); 
%reshape the raw data from STATA, SA-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
H1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, SA-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
H2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
P=[H1;H2];
%pick=find(P(:,2)==17 | P(:,2)==10 |P(:,2)==44 |P(:,2)==46 |P(:,2)==62 |P(:,2)==90);
%P=P(pick,:);
J1=dlmread(textFileName2na); %the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName2a); 
%reshape the raw data from STATA, AF-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
HH1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, AF-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
HH2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
PP=[HH1;HH2];
allsys=PP(:,end);
%Specify each page of FSTA and DDOM: fs and G
fs=zeros(size(P,1)*size(choice,2),T+5); 
G=zeros(size(P,1)*size(choice,2),T+4);
MS=zeros(size(P,1)*size(choice,2),T+4); %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation
eer=zeros(size(P,1)*ns,T+3);  %### UPDATED on 03/16/2017: include the simulated expected errors: ahaid, year, option given, prob@t=1, prob@t=2, ..., prob@t=100

for i=1:size(P,1) 
    tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
    exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
    exomkt=sortrows(exomkt,1);
    altemp=[exomkt;tempHH];  %group all the rival hospitals (standalone and affiliated) together (to do the random draw)
    tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
    alltemp=[tempmkt;tempHH]; %group all the hospitals (including the one considered) (standalone and affiliated) together (to do the random draw)
    %The following is to compute fs0 and g0. Then add [year, bdtot, fs0]
    %into fs and add [year,g0] into G.
    fs0=P(P(:,1)==P(i,1),4);
    mkts0=0;
    prechs=altemp(:,4);
    prechs(prechs==1)=[];
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
       mktsh=occur/sum(occur);
       mkts0=mktsh(fs0);
    end    
    prechs(prechs==2)=[];
    prechs(prechs==13)=[];
    if isempty(prechs)==1
       g0=0;
    end
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
     %  if occur(2)>0
     %     occur(2)=1;
     %  end
       occur=reshape(occur,1,[]);
       mktdom=choice(occur==max(occur));
       g0=max(fs0-mktdom==0);
    end
  %%%%%%%%%%%%%%%%%%%%%%%%%
  %consider the case where the hospital has already installed EMR  
   if P(i,4)>1   
        want=zeros(size(choice,2)-1,T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2)-1,T+1);  
        mkts=zeros(size(choice,2)-1,T+1);
        tempeer=zeros(ns-1,T+1);  %### UPDATED on 03/16/2017: a section of eer: the probability at each period
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market. ### 04/21/2020: should be
        %for stand-alone hospitals whose competitors are ALL stand-alone
        %hospitals, too.
        if isempty(sys)==1
          for l=2:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;
              if isempty(altemp)==1
                  newchs=l;
              end
              if isempty(altemp)==0
                 w=altemp(:,5:17);
                 cc=cumsum(w,2);
                 cc(:,end)=1;
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1  %for monopoly market
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l;
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                      %only for DEBUG
                      %newchs=[fs0;altemp(:,4)]; 
                      %newchs=[fs0;exouniv(:,4)];                     
                      %debug END                                          
                      
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                                        
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0;
                       %mkts(l-1,t+1)=0;
                     end
                     if isempty(prechs)==0
                        occur=histc(prechs,choice)';
                        occur=reshape(occur,1,[]);           
                        vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                        %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                        %vmshf=vmshff;
                        prechs(prechs==2)=[];
                        prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0; %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4);
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);  %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".                   
                     Tchoice=choice';                     
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);
                    
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)'; 
                     
                     % ### UPDATED on 07/15/2020: construct bdtot matrix
                     bdcol=tempmkt(:,3);                                                               
                     % ### UPDATED on 04/26/2020: generate yrsince
                     %yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:)...
                       +vmktf.*bdcol(:,ones(ns,1))*alpha(end);
                       %+vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:)...
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:)...                       
                          +vmktf.*bdcol(:,ones(ns,1))*tau(end);
                      % +vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:)...
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                                        
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);  %###!!!If the first one changed back to using max, remember to use "altemp"!
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];                   
                     %### UPDATED on 1/4/2017: re-define "selfnewchs" by
                     %figuring out what brings in largest payoff given the
                     %rivals' new choice, from "newchs(2:end)".            
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),2000+k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),(2000+k)*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market. ### 04/21/2020: consider hospitals whose competitors
        %include affiliated hospitals
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:); %find out the outside ones
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=2:13
              %w=altemp(:,5:17);
              w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1;
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l; 
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                      %only for DEBUG
                      %newchs=[fs0;altemp(:,4)]; 
                      %newchs=[fs0;exouniv(:,4)];                     
                      %debug END                                                               

                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     
                    alltemp(:,4)=newchs(1:size(alltemp,1),:);
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                     
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0; %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0;  %### UPDATED on 07/16/2018                     
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                  
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                    vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                    vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                    

                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);                           
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);                    
                     %bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);                     
                     %bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)>1);                                      
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                 
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                     % +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                     % +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                        
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';                    
                     % ### UPDATED on 07/15/2020: construct bdtot matrix
                     bdcol=univ(:,3);                     
                     % ### UPDATED on 04/26/2020: generate yrsince
                     %yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:)...
                       +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:)...                                         
                       +vmktf.*bdcol(:,ones(ns,1))*alpha(end);
                     %+vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:)...
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:)...                                                            
                         +vmktf.*bdcol(:,ones(ns,1))*tau(end);
                     %+vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:)...
                     v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:);
                      %+yrsince*Eta(6,:);                                                           
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:);
                      %+yrsince*Lambda(6,:);                                        
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                     
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);  %!!! change it to "univ" or "exouniv" depends on random or max.
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),2000+k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),(2000+k)*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
   end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %consider the case where the hospital has never installed EMR
   if P(i,4)==1  
        %tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
        %exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
        %exomkt=sortrows(exomkt,1);
        %altemp=[exomkt;tempHH];  %group all the hospitals (standalone and affiliated) together to do the random draw)
        %tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pc1 pc2...pc13
        %alltemp=[tempmkt;tempHH];
        want=zeros(size(choice,2),T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2),T+1); 
        mkts=zeros(size(choice,2),T+1);
        tempeer=zeros(ns,T+1); 
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market.
        if isempty(sys)==1
          for l=1:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;             
              if isempty(altemp)==1  %consider the case of monopoly
                  newchs=l;
              end
              if isempty(altemp)==0
                 w=altemp(:,5:17);
                 cc=cumsum(w,2);
                 cc(:,end)=1;                 
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                      %only for DEBUG
                      %newchs=[fs0;altemp(:,4)]; 
                      %newchs=[fs0;exouniv(:,4)];                     
                      %debug END                                                                
                      
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];
                                          
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;   %### UPDATED on 07/16/2018
                     end
                     
                     
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);     %### UPDATED on 07/16/2018                 
                     Tchoice=choice';      
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                   
                   curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);


                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);    
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);                  
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);                  
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';     
                     % ### UPDATED on 07/15/2020: construct bdtot matrix
                     bdcol=tempmkt(:,3);                                          
                     % ### UPDATED on 04/26/2020: generate yrsince
                     %yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:)...                       
                       +vmktf.*bdcol(:,ones(ns,1))*alpha(end);
                       %+vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:)...
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:)...                         
                         +vmktf.*bdcol(:,ones(ns,1))*tau(end);
                       %+vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:)...
                      %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;              
                     %selfnewchs=find(w(1,:)==max(w(1,:)));
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l,t+1)=newchs(1);  
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]; 
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];  %### UPDATED on 07/16/2018
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),(2000+k)*ones(size(tempeer,1),1),tempeer];
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market.
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that have belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:);
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=1:13
              w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1; 
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                      %only for DEBUG
                      %newchs=[fs0;altemp(:,4)]; 
                      %newchs=[fs0;exouniv(:,4)];                     
                      %debug END                                                                
                      
                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     alltemp(:,4)=newchs(1:size(alltemp,1),:);
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];                               
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;   %### UPDATED on 07/16/2018                   
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                 
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                     vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                     vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                   
                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);  
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);    
                                     
                    % bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);     
                    % bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)==1);                                        
                    % v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                    %  +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                   
                    % v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                    %  +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                    % v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                    %  +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                    % v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                    %  +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                    
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)'; 
                     % ### UPDATED on 07/15/2020: construct bdtot matrix
                     bdcol=univ(:,3);                     
                     % ### UPDATED on 04/26/2020: generate yrsince
                     %yrsince=ones(size(vmktf,1),1)*(tbar-t)*(t<=tbar)+zeros(size(vmktf,1),1)*(t>tbar);
                     %yrsince=ones(size(vmktf,1),1)*tbar;
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     % ### UPDATED on 04/26/2020: further include
                     % vdom*yrsince and yrsince*i.vid
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:)...
                       +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:)...                                         
                       +vmktf.*bdcol(:,ones(ns,1))*alpha(end);
                     %+vmktf.*yrsince(:,ones(ns,1))*alpha(9)+yrsince*Gamma(6,:)...
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:)...                                                            
                         +vmktf.*bdcol(:,ones(ns,1))*tau(end);
                     %+vmktf.*yrsince(:,ones(ns,1))*tau(9)+yrsince*Zeta(6,:)...
                     v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:);
                      %+yrsince*Eta(6,:);                                                           
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:);
                      %+yrsince*Lambda(6,:);
                                        
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;
                     %{
                     selfnewchs=find(w(1,:)==max(w(1,:)));
                     w(1,:)=[];
                     cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l,t+1)=newchs(1); 
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g];   
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),(2000+k)*ones(size(tempeer,1),1),tempeer];

        end
    end
end
G(:,[3,4])=G(:,[4,3]);  %exchange g0 with action_given
MS(:,[3,4])=MS(:,[4,3]); % ### UPDATED on 08/22/2018: exchange mkts0 with action_given
FSTA=[FSTA;fs];
DDOM=[DDOM;G];
MKTS=[MKTS;MS];
EERR=[EERR;eer];
end
dlmwrite(textFileName3,FSTA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName4,DDOM,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName5,EERR,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName6,MKTS,'delimiter', ',', 'precision', 8)


toc


